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ABSTRACT 



Ekman’s linear equations for time dependent flow (neglecting 
wind stress) are solved using a time dependent Green ! s function and 
the method suggested by Welander (1957). The solution represents 
the vertical velocity profile in terms of the local time history of 
the changes in sea surface elevation determined by the divergence 
of the flow in the vertically integrated continuity equation. 

A fully implicit finite-difference scheme is developed to re- 
present a time dependent seiche oscillating across a shallow infinite 
channel. The transients associated with the formation of the bottom 
spiral are clearly represented by the model and the influence of 
friction and Coriolis are individually and collectively introduced. 
The model allows independent calculation of velocity, volume trans- 
port, sea surface elevation, bottom stress, and the total energy 
balance of the system. The numerical scheme provides a method for 
adequately describing and investigating certain classes of time 
dependent motion, and its development suggests a mechanism for 
improving numerical prediction of storm surge. 
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I. INTRODUCTION 



"The purpose of computing is insight not numbers" 

(Hamming, 1962) 



A. BACKGROUND 

Oceanographers use a number of techniques for the study of 
oceanic flow including analysis of hydrographic data, use of 
hydraulic laboratory models, and mathematical treatment of the 
basic equations of motion. Initially oceanic circulation was 
treated as a steady state problem, and various investigators have 
been successful in representing and explaining many of the gross 
features of oceanic circulation using each of these techniques. 

Mathematical treatment of oceanic circulation began with 
Ekman T s analysis (1905) of wind-driven circulation and his results 
have been extended by a number of authors. (Sverdrup, 1947 > Munk, 
1950, Bryan, 1959) For the most part the study of large scale 
ocean dynamics has been done with steady state formulations. 

Mathematical treatment of time dependent motion is far less 
complete in the literature than that of steady state. It is also 
far more difficult to compare results due to the scarcity of 
reliable data. The complex nature of the solutions to time de- 
pendent boundary value problems further increases this difficulty. 
These solutions take the form of varying differential, integral, 
and infinite series equations whose numerical evaluation was not 
feasible until the advent of the high speed computer and the de- 
velopment of appropriate numerical techniques. 
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In the last decade these methods have been applied, through 
numerical models, to the problems of explaining and predicting the 
various details of fluid motion in both the atmosphere and the 
oceans. Although much success has been attained with these methods, 
many critical questions still need explanation. Many of these 
questions relate to the structure of the various frictional boundary 
layers, the details of their development, and the resulting effects 
on the interior flow. 

Many complicated engineering problems occur in near shore and 
other shallow water regions. The time dependent nature of these 
areas is well documented and the resulting flow is greatly modified 
due to interaction with the ocean bottom. Therefore an understanding 
of the nature of time dependent flow at all levels is fundamental to 
dealing with problems in these areas. The question of how readjust- 
ment takes place and how the many transients can be appropriately 
represented in the bottom frictional boundary layer has not been 
adequately explained. Further, it is necessary to understand the 
dynamic coupling between the flow and the development of the bottom 
boundary layer before accurate modeling techniques can be developed. 
It is towards this understanding that this work was directed. 

B. PREVIOUS RELATED RESEARCH 

Work by Welander (1957) dealt with obtaining time dependent 
solutions using Ekraan dynamics. His results were presented in 
terms of a time dependent Green’s function and were related to the 
problems of predicting tides and storm surges. He demonstrated 
that the velocity profile, volume transport, sea level elevation, 
and bottom stress could be exactly represented by his 
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"integr ©-differential" solution. It was suggested that this 
solution could be used to improve storm surge and tide predictions. 
A step-by-step numerical scheme was proposed, but the details of 
the scheme were not presented in the literature. 

Galt (1970) covered much of the same material that Welander 
presented but used a different method of solution. The solution 
allowed a more explicit formulation of results which provided more 
insight into the possible transients present and the relevant time 
scales of each. He discussed the need for information relating to 
time dependent flow in shallow water regions, and suggested that 
the numerical investigation of the interaction between a seiche 
and its associated transient bottom spiral might be useful in de- 
termining how the various transients present are effectively 
coupled. 

C. THESIS OBJECTIVES 

Since the major area of interest was to be the bottom boundary 
layer, the first objective was to formally solve the governing 
equations neglecting wind stress. The solution would therefore not 
allow a surface Ekraan spiral to develop and effect the interior 
flow. It will be showed that the solution can be represented by 
two equations, one integral and one differential. These equations 
are similar in form to the ones developed by Welander and Galt. 
Although both of these authors proposed that these equations could 
be solved numerically using step-by-step procedures, it was not 
initially apparent whether a stable and accurate numerical scheme 
was feasible. Therefore the second objective was to formulate a 
numerical model capable of simultaneously solving these equations. 
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Finally, as suggested by Galt, this model would be applied to 
represent a seiche to examine the relationships which developed 
between the various transients. 
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II . THE FORMAL SOLUTION 



A. GOVERNING EQUATIONS AND ASSUMPTIONS 

The momentum equations considered are essentially those used 
by Ekman. The non-linear and horizontal friction terms are assumed 
small and neglected. The Coriolis parameter, the density, and the 
coefficient of vertical friction are all assumed constant. Thus: 



bu 

bt 



fv 




bv 

bt 



+ 



fu 




g 



b£ 

bx 
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$1 

by * 



( 1 ) 

( 2 ) 



Where u and v are velocities in the x and y-dir ect ions , f the 
Coriolis parameter, g the acceleration of gravity, and k the 
constant coefficient of vertical friction. §(x,y,t) is the ele- 
vation of the free surface and the z-axis is considered positive 
up . 

The third equation necessary to solve for u, v, and 5 is the 
vertically integrated continuity equation. 



b| bu bv 

bt bx by ~ 



(3) 



where 



r 0 

U = \ u dz 
J -h 



c 0 

; V = \ v dz 

-h 



( 4 ) 



is the volume transport and h = h(x,y) is the depth of the water. 
§ is assumed small compared to h. 
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Neglecting the effects of wind stress on the motion results in 



the following linearized boundary conditions on u and v: 



(fe) 2 =o = ° ’ (& 2 =o = 0 



(u) _ = 0 ; (v) _ = 0 . 

V ' X=-h v * = -h 



z = -h 



( 5 ) 

( 6 ) 



For this development it was assumed that initially a state 
of rest exists: 

(u) t=0 = 0 ; (v) t=0 = o . (7) 

As in Welander T s work the following complex notation was 
introduced: 

w = u + iv 



bn bx by 

Where i = V - 1 . Using this notation eqs . (1), (2), (6), and (7) 
can be expressed as follows: 

bw . b 2 w . , b| 

Bt ■ k rr ^ = • s s 

oz 

f— ^ = o 

( b 2 ; 2 =o 



( 8 ) 



< w) t=o = 0 

< w ) 2 =-h = 0 



( 9 ) 

( 10 ) 
( 11 ) 
( 12 ) 



Equations ( 9 ) through (12) define the initial/boundary valued 

v e 

problem that must be solved to specify u and v. Although <5S> in 
eq. (9) is not externally specified it can theoretically be deter- 
mined by integrating eq. ( 3 ). This makes it possible to treat it 
as a time dependent forcing term (i.e. ^ (t)). 
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B. METHOD OF SOLUTION 



To solve the problem begin by multiplying eqs . (9) through 

if t * if f 

(12) by e and introducing the relation w = we • Then eqs. 

(9) through (12) can be written as: 



bw . b w ift b§ ,io\ 

k — . . ge _ (t) (13) 
02 


<5f\=o - 0 


(14) 


* 




o 

II 

o 

II 


(15) 


(w*) _ = o . 

' ' z=-h 


( 16 ) 



A solution to eqs. (13) through (l6) can be obtained by using 
a time dependent Green’s function approach as described below. 
Assume 

W * = 0 m (t > Z m (z ) (17) 

where Z (z) satisfies the boundary conditions and 0 (t) satisfies 
the initial conditions. Z m (z) is the eigen-function determined by 
applying standard separation of variables to the homogeneous part 
of eq. (13) using eqs. (14) and (l6). It is found to be: 





Z m (z) = cos(pz) 


( 18 ) 


where 




D = ( 2m -i) 

F 2 h 


(19) 



0 ra (t) is determined by substituting eq. (17) into eq. ( 13 ) 



requiring that 0 m (t) satisfy eq. ( 15 ), and using the orthogonality 



of the cosine function to determine the appropriate constants. 
Thus 

a (t\ = — 

ph 






( 20 ) 
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Therefore a solution to eqs . (13) through ( l6 ) is 

* , , , 2 g ^ ( -1 ) m cos(pz ) ~ k P 2 ( t-t ' ) -if t’ b| , ^ x 

w (x,y,z,t) = -jS E =1 J > 'P' \ e p >e (21) 



For the general time dependent solution of eq. (9), dividing 
eq. (21) by e 1 ^* yields: 



-(x.y.x.t, ■ g L ^^- agl S t *- <l,>2 * lf)<t - t,) ^ (t ' )dt ' (22 > 



C. DISCUSSION OF THE FORMAL SOLUTION 



It is apparent at a glance that eq. (22) satisfies the boundary 
and initial conditions. Substituting eq. (22) into eq. (9) yields 
as a residue: 



S=i 



Izll 



m 



ph 



cos(pz) = 



(23) 



It is easily verified that the left hand side of eq. (23) is 
simply the cosine expansion of the right hand side, and therefore 
eq. (22), in the limit, is an exact solution to the initial 
problem. 

Equation (22) represents the horizontal currents at all depths 
resulting from the pressure gradient produced by the slope of the 
sea surface at a particular point and instant in time. The integral 
term carries a running summation of the currents resulting from the 
slope in the past history of the flow. These currents are modi- 
fied in time by a rotation and decay represented by the weighting 
function 

e -(kp 2 + if)(t-f). 

This weighting function introduces two interesting time scales 
into the problem. The first is the inertial rotation period, which 



15 



varies as a function of latitude. For raid-latitude locations it is 
on the order of 20 hours. The second time scale is introduced by 
the exponential decay (T^) represented by 

T = _1_ = 4h 2 

d kp 2 k(2m-l ) 2 tt 2 

This is a function of depth and friction, but for typical values 
2 

(k = .01 ra /sec, h = 50 ra) it is on the order of 2 8 hours for the 
first terra in the series to decay to 1/e of its initial value. 

This is the slowest terra to decay. 

The motion represented by (22) depends on the surface slope 
and thus new time scales are introduced into the solution. These 
depend on the horizontal dimensions represented and enter through 
the integration of the continuity equation (3). Time scales re- 
lated to set up of the sea surface and seiche periods will govern 
the time rate of change of the surface slope. This change 
initiates the transients present in the geostrophic and bottom 
currents, which show the inertial rotation and exponential decay 
previously noted. 
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III. THE NUMERICAL SOLUTION 



A. FORMULATION OF THE MODEL 

1 . Establishment of a Representative Cross Section 

In formulating the numerical scheme it was decided to 
simplify the model by neglecting the derivatives along the y-axis. 
To accomplish this without the loss of any essential physics, it 
was decided to allow the seiche to oscillate across a narrow, 
shallow channel of infinite length. This prevents sea surface 
slopes from developing along the length of the channel but is 
otherwise non-restr ictive. The simplification further reduces the 
problem to one of looking at a cross-sectional distribution of 
variables in the channel. 



a cross-section of constant rectangular dimensions. The width and 
dqpth of the basin considered were fixed at 5000 meters and 30 
meters respectively, and Figure 1 is a pictorial depiction of this 
cross section. 

These simplifications allow eqs . (3) and (22) to be re- 

presented in the model as 



It was further decided to restrict the model by considering 




(25) 




-(kp 2 +if)(t-f) b£ 
bx 



t ' ) dt ' ( 26 ) 



17 



30 


m 


BASIN CROSS-SECTION 






1 5000 m 1 




Z:o 



z:-h 

XrO X=L 



MODEL CROSS-SECTION 



Figure 1. Pictorial Depiction of Model Basin 
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2 . Establishment of Governing Non-dimensional Ratios 



To represent the various time scales and frictional co- 
efficients possible with more general time dependent motion the 
model was scaled to the period (T) of a free oscillating seiche in 
the model basin, and the coefficient of vertical friction determined 
by Ekraan's "depth of frictional resistance (D)". Equation (26) was 
scaled as follows: 



T 



2L 

(gh 



D 




(27) 

( 28 ) 



Let, 



* 

t = Tt 

* 

t' = Tt* 
dt' = T dt'* 



Lx 



x = rrr (eleven grid points across the channel) 



&x 



10 

10a 



o 

■» 

bx 



* * * * 

where t , t' , | , x , are non-dimensional quantities and A is 

o 

the free surface elevation at x = 0, t = 0. Substituting eqs. (27) 
through ( 29 ) into eq. ( 26 ) yields: 



w 



= 2 a “ (— l) m cos(p 2 ) 10A o T + if ) ( t-t 1 ) *T 

h m=l p L d Q 2tt 



X »L 

bx* 



dt ’ 



(30) 



Recalling that p = — — - — , and O' = — eq. (30) may be written 
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40g 



W = 



A T 
o 



E -i 

ra -1 



124*. 

(2ra-l )n ^0 



. .m 

•1) cos 



( 210 - 1 )' 

*r 



(-) 



(t-t 1 ) 




bx 



( 31 ) 



It was determined that the two non dimensional ratios, 
(f/CT) and (D/h), appearing in eq. (30) govern the time/spacial 
scale and frictional dependence of the motion. For discussion and 
computation the first was designated (PCF) and the second (PCD). 

The first of these ratios represents the relationship 
between the earth T s inertial frequency and the frequency of the 
basin in the following manner. 

f = 2 O s in 0 

where Cl is the frequency of the earth's rotation and 0 is the 
latitude . 

O’ = 2 rr F , 
ra 



where F is the resonance frequency of the basin. Since no external 
ra 

forcing terms were used in the solution, F^ is also the approximate 
frequency of the motion. Thus: 



— = ^ sin0 
O “ F 

m 

The model was scaled for (0 £ PCF ^ 1) . A value of 
(PCF = 0) represents motion of a time scale which occurs too 
quickly to be affected by Coriolis accelerations, while a value 
of (PCF =1) represents motion whose time scale is equal to that of 
inertial rotation and subsequently greatly influenced. 
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The second term, PCD, is the ratio of the depth of 
frictional influence (Sverdrup, Johnson, and Fleming, 1942) to the 
actual water depth. The coefficient of vertical friction was de- 



termined in the model from eq. ( 28 ) as 



k 



(PCD X h) 2 X f 



( 31a) 



The model was scaled for (0 ^ PCD ^2). A value of (PCD = 0) 



represents a very deep basin where frictional effects on the motion 
are minimal, while (PCD = 2) represents a basin where friction domi- 
nates the motion. 



model is controlled strictly by the resonance of the basin, then 
the spacial scale of the basin is determined by a choice of PCD 
and PCF . PCD will govern the approximate depth of the basin while 
PCF determines the basin width through eq. (27). Certain con- 
clusions about motion caused by external forces (i.e. tides, 
atmospheric pressure) producing time scales independent of the 
basin size were also inferred from the model and are discussed in 
the conclusions. 

3* Adaptation of the Formal Solution to the Model 

To establish two equations that could be solved simulta- 
neously for sea surface elevation and transport, eq. ( 26 ) was 
integrated in the vertical using eqs . (4) and (8). Thus: 



and W - U + iV. Equation (32) represents the total volume trans- 
port associated with the flow at a particular point in the channel. 



If the time scale of the motion being represented by the 




+if)(t-t') 

£>x 






(32) 
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The real part of eq. (32) represents the cross-channel transport, 
while the imaginary part the along-channel transport. 

The numerical scheme was established to solve eqs . (25) 
and (32) for § and W. The velocity profile, bottom stress, and 
energy balance associated with the flow were also recoverable from 
the scheme, and the details for this recovery will be discussed 
later . 

4 . Initial Conditions 

The model assumes an initial state of rest in the motion 



and an 


initial set-up of 


the free 


surface . 


This was represented 


as 


W(x,0) 


= 0 , 




(33) 


and 


5(X,0) 


= A cos 
o 


(TTX) 

L ’ 


(34) 


wher e 


W = W(x,t) and | = 


?(*, t) . 







5. Boundary Con dit ions 

Since the governing equations were solved without reference 



to horizontal boundaries, these conditions had to be built into the 
model. The model therefore assumes no flow across the horizontal 
boundaries, or 

W(0,t) = W(L, t) = 0. (35) 

These conditions were satisfied by requiring the surface slope to 
remain zero at these boundaries at all times. Therefore eq. (32) 
is forced to meet the conditions of eq. (35). 

6 . The Step-by-step Scheme 

Both Welander and Galt suggested that a step-by-step 
numerical scheme could be used to solve eqs. (25) and (32). As a 
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first solution the staggered central difference scheme suggested 
by Galt was used and is summarized below. 



t 



At Given initial conditions on W. 



t = 0 Given initial conditions on 

t = At Calculate W from eq. (32) using ^ at t = 0. 

t = 2At Calculate 5 using an explicit finite-difference 



scheme on eq. (25) and U at t = At. Calculate a new 



surface slope. 



t = 3 At 



. . etc. 



The term At is the time step of the explicit scheme. 

This scheme was developed neglecting Coriolis (f = 0) in 
eq. (32). The scheme proved to be stable and accurate subject to 
the following stability criterion: 



The expression Ax is the spacial step across the channel. It was 
felt however that a more precise solution, not subject to such a 
restricted stability criterion, could be formulated using a totally 
implicit finite-difference scheme. The development of this method 
is now presented. 

7 • The Implicit Finite-difference Scheme 

The implicit scheme used is a modification of a method 
described by Richtmyer and Morton (1967). 

Start by approximating the time rate of change of the 
transport as follows: 




(36) 



bw _ < W W " < W >t 



(37) 



bt At 
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From eq. (3 2 ) 



&W _ _I_ (~2£ v _L C t + At e " (kp2 + if )( t+At ~ t ' ) $if t 'ldf 

bt ■ At l h a=i 2 [i 0 e bx (t )dt 



_ e -<xp 2 *if)<t-f> || ( t.)dt 



( 38 ) 



-t pt+At 

By expanding the first integral as \ + \ and regrouping terms, 

J 0 J t 



eq. ( 38 ) may be written as 






i e -(T 2 *«)(t*«-.') 5£.(t*)df 

At ^ h m=l 2 f)x v 7 



(39) 



Eqs . (25) and (32) can then be represented in the following form: 



( £L) = _ I 

'bt' A_t 2 
2 



bU\ / bu 

^t ^ t+At 



(40) 



= ( ci * + ,c£\ 

bt + At At ' ' 2 At 

T 2 



s), * (s) 



t + At 



(41) 



wher e 



2 . 



C1 * = ~~h~ m=l S " (kP +if)(t " t,} (e" (kp +if ) At -l) §|(f)df ( 42 ) 

P 0 



«* . -22 5 1 

h m-1 2 

P 



f>t+At e “ ( kp 2 +if ) (t+At-t 1 ) dt f 
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Equations (40) and (4l) can be represented in Matrix form as 
follows : 



( 5 T> ttA t =A *f^ [< 6 >t * < e W] 



( 43 ) 



wher e , 




\ 



A = 



{ Cl/At 

T ! 



B = 



0 0 C2/At 

0 0 D2/At 

-10 0 



and. 



Cl = Real { Cl*} ; C2 = Real {C2*} 
D1 = Iraag {cl } ; D2 = Iraag { C2 } 



(Real { } and Iraag { } represent the real and imaginary part of 
the argument in braces) . 

The following implicit finite-difference equation was used 
to represent eq. ( 43 ): 



9 n+1 - 9 n 

_J 2 . 

At 



= A + 



B 

4a> 



- er} + e" - e n , (W) 

3+1 3 -I j+l 3-1 ' 



where 9 . 

J 



n _ qU At 



3 Ax 
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and n = (1,2,3, . maximum number of time steps) 






are the grid points in the time and x-directions respectively. 
Letting 



O' 



AtB 

4ax 



and 



d. = ao" . +©" - . + AtA, 

J D + 1 J 3-1 



eq. (44) becomes 



- oteV} * e n+1 * as 1 ?*; 

j+i j j-i 



D. . 

D 



(45) 



The left hand side of eq. (45) contains the unknown values of 
(9) to be calculated, whereas the right hand side ( D_. ) are all 
known values . 

A recursion relation of the following form was assumed. 



e n+1 = e. o n i; * F. 

J 3 3+1 3 



where 



and 



E. 

J 



F . 

J 



e ll 


e 12 


e 13 


4 


4 


e 23 


4 


4 




( 4 ) 






) 








(46) 
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If eq. (46) holds for all 9 ^ , then the model's left hand boundary 
conditions require that 



’ 0 0 0 


f° N 


0 0 0 


and F = \0 1 




° 


0 0 1^ 


i»J 



(47) 



Substituting the expression for 
eq. (45) yields: 



9 n+1 

j-l 



derived from eq. 



(46) into 



9 n + i 

J 



a 



I4 * E j-iJ 



9 n+ i + 

j + 1 



D.-a f. , 
2 J" 1 

i+a. e. , 

J- 1 



(48) 



Equating the right hand side of eq. (48) to eq. (46) gives the 
following recursion relationship: 



E . 
J 



F. 

J 



CL 



i+a e . . 
j-i 

D.-O! F. ‘ 

J 111 



I-KX E 



3-1 



3 s i; 
3 * 1; 



(49) 

(50) 



and I is the ( 3X3 ) identity matrix. 



Using eqs . (47), (49), and (50), values of E^ and F ^ were 
calculated inductively in order of increasing j(j = 1,2,3,..., J-l). 
Since the value of 9^*^ is given for j = J-l by the model's right 



boundary condition (U*J = = 0 and = §*! ), values of 

<J " J- 



hand 

n+ i 

9_. were then calculated inductively from eq. (46) in order of 



decreasing j ( j = J-l, J-2, ... , 1). This completes the calculation 

.. , J-l). The details of the mathematics for 



of 9*? +1 (j = 1,2,3, 



representing E^. , F^ , D^. , and eq. (46) are presented in appendix A. 
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* * 

8. Numerical Treatment of Cl and 02 

* * 

The numerical scheme for representing Cl and C2 is 

extremely lengthy and is merely outlined in this section. For a 

more detailed presentation the reader is referred to appendix B. 

* 

The integration of Cl was performed by dividing the total 
integration time into small time intervals (At) and assuming that 
the surface slope and effects of inertial rotation could be re- 
presented adequately by a constant throughout j:he interval. These 
values were calculated at the mid-point of each At interval and 
used in the scheme as representative values throughout the entire 
sub-interval. This permitted the integral to be evaluated analyt- 
ically over each interval with a maximum error on the order of 

* 

these approximations. Using these assumptions Cl was numerically 
represented as 



Cl 



_ MAXM, 

= "2a s / JL 

h j 2 


^(l-e- k P 2at )' 


i-l r 

S 


M=1 / p 


.kp 


N=1 • 



‘bx L (e 



kp at -if At , . 

e -1) 



(e -if(I-VN)At )(e -kp 2 (I-l-N)At ) 



(51) 



for I = 2, 3j ••• y and Cl = 0. The value of the integer I ranges 
from 1 to MAXT (maximum problem time) and iAt represents the elapsed 
time from initial conditions. L takes on integer values from 1 to 
MAXL (maximum number of cross-channel grid points) and LAx is the 
distance from the left boundary. The integer M ranges from 1 to 
MAXM, and determines the maximum number of modes in the series that 
are included in the approximation of Cl • It was determined from 
eq. (23) that MAXM = 25 would represent the solution with a maximum 
error of three percent, and this value was used throughout all 



28 



numerical evaluations. The integer N takes on the values from 1 

to (1-1) and determines, through the exponential terms, how much 

rotation and decay has taken place in the recent time history of 

each transient considered. Since the terras of the series decay 

more rapidly with increasing friction (especially the higher modes), 

it was not necessary to carry their integration as far back as 

N = 1. It was decided to include only the terms of significant 

magnitude, and this was accomplished by requiring that when the 
2 

value of kp (I-l-N)At in eq. (51) was greater than 10, the cor- 
responding terms were discarded. 

* 

The integration of C2 involves only one time interval 

(t to t+At). Again the amount of inertial rotation was evaluated 

at the mid-point and assumed constant throughout the interval, and 
* 

C 2 was evaluated as the following constant: 



C2 



* = “2fl r™ e " if(?sAt ) 

h M=1 



■±5 (l-e- k P 2At ) 



. k P 



(52) 



* * 

Cl and C2 both contain the terra 



-±rr (l-e' kp il ) 



kp 



, and for k = 0 this terra can 



not be numerically evaluated in this form. However this difficulty 

2 

was overcome by expanding e ^ in a Taylor series, i.e. 



(l . kp 2 st * Us!«l 2 . 



( 53 ) 



and substituting this series for the term in eqs . (51) and ( 52 ) 

whenever the value of kp 2 At was less than 10 - ^. 

The method for separating the real and imaginary parts of 
* * 

Cl and C2 involves using the relation 

e = cos(ft) - i sin(ft) 
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in eqs . (51) and (52) and gathering terras. The details are 
covered in appendix B. 



9 . Velocity Calculations 

The velocity profiles across the channel were calculated 

by the model at discrete time intervals determined by an input 

parameter. Noting the simularity between the integral of eq. (26) 

* 

and Cl , it was numerically possible to compute and store the terms 

of these integrals needed for velocity, during the calculation of 
* 

Cl . For velocity profile computation, these stored values were 
recalled and applied in eq. (26). The details for numerically re- 
presenting eq. (26) is also presented in appendix B. 

10 . Bottom Stress Calculations 

The bottom stress (T) associated with the motion was 
calculated using the following relation: 



of eq. ( 26 ), bottom stress is easily calculated whenever the 
velocity profile is recovered. Appendix B also contains the 
details for this calculation. 

11 . Energy Calculat ions , 

The total energy balance associated with the seiche is a 
sura of its potential and kinetic energy. This balance was 




or fr om eq. ( 26 ) , 




-(kp 2 +if)(t-t') b£ 
bx 






( 54 ) 



Since the integral term of eq. ( 54 ) is identical to that 
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calculated as follows: 




( 55 ) 

( 56 ) 



PE and KE are the potential and kinetic energies and (p ) is the 
is the density of the water. Equations ( 55 ) and ( 56 ) were evaluated 
using Simpson f s 1/3 rule (McCalla, 1967) whenever velocity was 
calculated . 

The energy calculations were monitored during model 
testing as a calculat ional check since small numerical errors 
were quickly apparent in increasing energy levels. 

12 . Scaling of Model Output 

Since the governing equations for the model were scaled 
to produce the ratios PCF and PCD, the output of the model must 
be appropriately re-scaled to obtain real values for the variables. 
The series of graphs were produced on a CALCOMP plotter with cor- 
responding values from the model. The various axes were auto- 
matically scaled by the computer to fit the data being represented. 
These axes were then re-scaled using a non-dimensional scaling 
factor derived from the continuity equation to allow real values 
to be easily obtained. The exception is the figures depicting the 
energy balance, which are not scaled and should only be used to 
indicate relative changes in energy levels. 
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The scales for the various axes are as follows: 



AXIS 



SCALE 



Surface Elevation 



Time 



Velocity 



%/k 

o 

h 

t (gh) 2 

2b 

2w (h/q)^ 



Transport 



Bottom Stress 



2W 



C 2 A o (gh)' 



2T(h) 



3/2 



C 3 A o kg3 



The values of g, h, k, A q , L are determined from the 
actual basin being represented and are scaling con- 

stants dependent on the parameters of the prototype basin and the 
scaling of the graphs by the computer. These values are calculated 
for each figure and inserted in the scale. The variables, t, w, W, 
5, and T are the real values associated with the basin being re- 
presented and may easily be obtained from the graphs by entering 
appropriate arguments for g, k, h, A q , and L in the corresponding 
scale. 



B. TESTING THE MODEL 

1 . Establishing a Reference Solution 

One of the major difficulties in numerical analysis is 
verifying the accuracy and stability of the scheme. Since many 
numerical solutions deal with problems that have no analytical 
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solutions it is common practice to test the accuracy of these 
models against empirical data, or analytical solutions involving 
degenerate cases. As previously noted, data pertaining to time 
dependent flow is extremely scarce, and it was decided to test the 
accuracy of the model against a known analytical solution for the 
degenerate case where PCD = 0 (i.e. no friction). Analytical 
solutions to this case are readily available in the literature and 
a modification of one proposed by Prodman (1952) was used in the 
verification of the accuracy of the numerical scheme. 

2 . Accuracy and Stability of the Model 

The accuracy of the model was determined by comparing the 
numerical representation of the free surface elevation at (x = 0), 
the raid-channel transport, and the total energy balance, with the 
analytical values. Figures 2 through 5 highlite this comparison 
and indicate that the accuracy of the model was a function of time 
step (At) and the Coriolis parameter (f). It was noted that as the 
time step was reduced the numerical solution appeared to converge 
to the analytical solution. The phase shift in surface elevation 
over five periods was reduced from 45 degrees for (At = T/20) to 
24 degrees for (At = T/40) and appeared to be independent of the 
value of Coriolis. This error might possibly be introduced into 
the system by approximating the sea surface slope as a constant 
value in the time history integration scheme. As (At) increases 
this linear approximation is extended over a greater time span 
and the resulting errors amplified. As the value of the Coriolis 
parameter was increased, the error in representing the amplitude 
of the free surface and total energy balance also increased, as 
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Analytic 

Model (At = T/20) 

Figure 2. Comparison of Free Surface Elevation for Model and Analytic Solution 
(PCF = 0, PCD = 0) . 
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Figure 3. Comparison of Free Surface Elevation for Model and 
Analytic Solution (PCF = 1, PCD = 0). 
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Figure 4 . Comparison of Energy Balance for Model and Analytic Solution (PCF=0, PCD=0) . 



o o 
cm <r 




4 -> 







04 CM 



Q) 

B 

•H 

H 

TJ 

Q) 

N 

•H 



CX5 

B 

u 

O 

2 



* 37 



Figure 5 - Changes in Energy Balance as a Function of Time Step 



depicted in Figures 3 and 5- The amplitude deviation was a maximum 
of one percentage for (At = T/20) and the energy rise was 8.5 and 
3.7 percentage for (At = T/20) and At = T/40) respectively. Errors 
of this nature did not appear for solutions which neglected Coriolis, 
and possibly were introduced by approximating the amount of rotation 
[i.e. e ^ ] as a constant in the integration scheme. It is 

not exactly clear how this error increases the energy in the model 
but it was observed that it appeared linear and very systematic. 



numerical instability, a standard Fourier series method (Smith, 
1965) of stability analysis was attempted on eq. (44)- However, 
due to the complicated nature of the amplification matrix, the 
eigenvalues could not be readily obtained and the analysis was 
abandoned. As an alternate test, the model was run with (At = T/5) 



further linear rise in the energy level. These errors were there- 
fore judged to constitute a problem in numerical accuracy and not 
stability, for the range in which the test cases were run. 



could be predicted and controlled by varying the time step, and the 
model was run maintaining a maximum allowable error of one percent 
in the elevation and transport amplitude, and four percent rise in 
the total energy balance. This was accomplished by using a time 
step of (At = T/20) for (PCF < .5), and (At = T/40) for (PCF ^ . 5 ). 
Since the errors in phase shift were strictly a function of time 
step, they were completely predictable and accounted for during 
data analysis. 



To determine whether this energy rise constituted a 



for (t = 5T). This corresponds to 




, and resulted in a 



It was further found that the magnitudes of these errors 
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It should be noted that for cases where friction is 



neglected (k = 0), the transients do not decay with time and must 
be carried throughout the integration scheme. Additionally, the 
bottom boundary condition degenerates to a free slip situation and 
the approximation of a constant velocity profile with a truncated 
cosine series becomes more inaccurate. The number of calculations 
increase rapidly as friction is decreased, which may introduce 
increasing errors in round-off. 

Based on the comparison with the analytical solution it 
was concluded that the accuracy of the model could be sufficiently 
controlled such that the details of time dependent motion, as 
described by the governing integral and differential equations, 
could be adequately represented and examined. 

C. MODEL LIMITATIONS 

The most serious limitation imposed on the model were the 
basic assumptions of Ekman dynamics used in the formal solution 
of the governing equations. These assumptions were considered 
necessary to simplify the mathematical treatment of these equations 
to the extent that a numerical solution was feasible. The additioanl 
simplifications imposed during the formulation of the model further 
limit its capabilities. It was felt however that valuable experi- 
ence in the techniques of numerically representing time dependent 
motion would be gained during the formulation of this simplified 
model, and insight gained from the numerical product. With this 
information available it might be possible to further refine the 
model and remove many of the more important restrictions. The 
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I 

limitations inherent in the present model are listed in this 
section and the more important ones will be discussed in section 
(V) with recommendations for improving the scheme. 

The basic assumptions made in arriving at eq. (22) have the 
following limiting effects of the model’s capabilities: 1) the 

non-linear terms in the motion have been ignored; 2) a variation 
in the coefficient of vertical friction in the layer was not 
permitted; 3) the assumption of constant density neglects the 
effects of stratification on the pressure gradients in the layer, 
and the subsequent modification of the vertical momentum transfer 
and 4) horizontal friction was considered insignificant and 
neglected. 

Simplification of the model resulted in the additional 
limitation: 5) the assumption of a constant basin neglects the 

important effects of bottom topography on the flow. This is es- 
pecially important for considering motion in shallow near shore 
areas where these effects are large. 
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IV. PRESENTATION OF RESULTS 



A. TEST CASES AND GEOPHYSICAL SIMULATIONS 

1 . Description of Test Cases 

In order to gain maximum insight into the physical processes 
occurring in the motion, as well as to study the effects of inertial 
rotation and friction on the variables, the model was run for 25 
test cases in the ranges of (0 ^ PCF ^ 1) and (0 ^ PCD ^ 2). 

Table I summarized the various controlling parameters of each run 
and includes the Central Process Unit (CPU) time used on the IBM 
360/67 system for each program. 

2. Description of Geophysical Simulations 

As interesting dynamic interactions developed in the test 
cases, it was decided to simulate actual geophysical erabayments 
to determine to what extent these phenomena were present and 
significant. The actual length, depth, and latitude of the basins 
were measured and used to calculate (f), (O’)? and (PCF). An Ekman 

depth of 20 meters was assumed in the calculation of (PCD) . Runs 
26-30 represent these simulations and the various imput parameters 
are summarized in Table II. 

B. DISCUSSION OF TEST RUN RESULTS 

1 . Frictional Effects on the Flow 

To describe the conditions where friction totally dominates 
the motion, the model was run with varying amounts of friction and 
(PCF =0). It must be noted that by neglecting the effects of 
inertial rotation, the Ekman depth is no longer specified as 



Table I. Parameters of Test Cases (runs 1-25) 



RUN NO. 


PCF 


PCD 


CPU( MIN-SEC) 


1 


.25 


.5 


6-2 


2 


.25 


1.0 


3-51 


3 


.25 


1.5 


3-8 


4 


.25 


2.0 


2-49 


5 


0 


.5 _ 


5-11 


6 


0 


1.0 


3-20 


7 


0 


1.5 


2-43 


8 


0 


2 . 0 


2-24 


9 


0 


0 


13-33 


10 


.25 


0 


14-58 


11 


.5 


2.0 


6-25 


L2 


.75 


2.0 


6-0 


13 


1.0 


2.0 


5-48 


14 


.5 


1.5 


7-25 


15 


.75 


1.5 


6-42 


16 


1.0 


1.5 


6-25 


17 


.5 


1.0 


9-40 


18 


.75 


1.0 


8-21 


19 


1.0 


1.0 


7-4o 


20 


,5 


.5 


16-14 


21 


.75 


.5 


13-52 


22 


1.0 


.5 


12-29 


23 


.5 


0 


62-48 


2 4 


.75 


0 


62-38 


25 


1.0 


0 


62-51 



NOTES: For the model basin 
10 T = 583.2 Sec 

2) CT = . 108 X 10 _1 Sec 



42 



Table II. Parameters of Geophysical Simulation 
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NOTE 



evidenced by eq. (28). Therefore, the frictional coefficients 
generated by runs 1 through 4 were artificially imposed on the 
model with (PCF = 0), and resulted in runs 5 through 8. 

The results for these cases are depicted in Figures 6 
through 8 and are simply a seiche damped by friction. As the 
friction coefficient increased the damping of the motion increased 
with a resulting increase in the period of oscillation. For the 
maximum case of (PCD = 2), the seiche was completely damped after 
4^ periods. These results were completely predictable and are 
presented to show that the model's representation of motion which 
is dominated by friction conforms to previous observations. 

(Pierson and Neuman, 1966) 

An interesting coupling between bottom stress and transport 
developed during these cases and a typical example is depicted in 
Figure 9. In all cases examined a phase lag quickly developed be- 
tween the transport and the bottom stress. For the case showed, 
the value of this lag averaged 45 degrees through the periods of 
oscillation. This phase lag can be explained by recalling that 
the transport is an integration of the motion over the entire water 
column, whereas the bottom stress is solely dependent on the 
velocity gradient at the bottom. Since friction has its greatest 
effect at the bottom, changes in magnitude and direction of the 
flow nearly always initially occur next to this boundary and prop- 
agate upwards in the column. The result is that changes in trans- 
port always lag changes in the bottom velocity gradient and 
consequently the bottom stress. This phase lag has obvious effects 
on the accuracy of commonly used numerical schemes where bottom 
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Figure 6. Effects of Friction on the Free Surface (PCF = 0) . 
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Figure 7. Effects of Friction on Mid-channel Transport (PCF 
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Figure 8. Effects of Friction on Energy Balance (FCF 
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stress is parameterized in terms of the transport. This technique 
places strong restrictions on possible interactions between the 
flow and the bottom, and although it may be sufficient for approx- 
imating steady state problems, it is unlikely that it will represent 
transient motion very well. 

Figure 10 represents typical velocity profiles for motion 
which is totally dominated by the effects of vertical friction. 

The increased period of oscillation and the mechanism causing the 
phase lag between the transport and the bottom stress is clearly 
evident. At times (t = T/2 ) and (t = 3T/4) the profiles are in 
transition from being totally positive to totally negative, and the 
transport and bottom stress are 180 degrees out of phase. 

These velocity profiles are typical of those found in other 
boundary layer flow and are similar to the wind profiles observed 
at the air sea interface over smooth water. They indicate that the 
stability of the flow has a strong frictional dependence. At 
(t = T/4) and (t = T) the profiles are maximum and, except initially 
in the upper layers, friction exerts a strong influence throughout 
the column. The result is a smooth regular flow and a well de- 
veloped frictional boundary layer with the region of greatest shear 
at the bottom. During the transition phases (t = T/2 to t = 3T/4) 
the profile changes direction and the water column appears to be 
less stable. The velocities in the column are minimum just prior 
to and after reversal and thus the influence of friction in the 
column is reduced. This region of reduced frictional influence 
moves up the column at the point where the flow changes direction. 

As the flow completes its reversal and the velocities in the column 
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Figure 10. Effects of Friction on Mid-channel Velocity 
(PCF = 0, PCD = 2) . 
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Normalized Depth (-z/h) 



increase, the effect of strong frictional influences is quickly 
re-established throughout the layer, the region of maximum shear 
shifts to the bottom, and smooth flow once again prevails at time 
(t = T). 

2 . Coriolis Effects on the Flow 

To gain insight into how inertial rotation modifies the 
transient motion, runs 10, and 23 through 25 were made neglecting 
friction (PCD = 0). 

It is important to realize that the relationship between 
the transport and the free surface elevation for a non rotational, 
frictionless seiche, as seen in Figures 6 and 7, is similar to the 
velocity/elevation relationship of the frictionless pendulum. At 
(t = 0), (t = T/2) , and (t = T) , the energy of the seiche is en- 
tirely potential, whereas at (t = T/4) and (t = 3T/4) it is 
entirely kinetic. At intermediate times the energy balance is a 
combination of potential and kinetic and is totally conserved. The 
free surface elevation and transport periodically oscillate between 
maximum positive and negative values of respectively equal magni- 
tudes, always maintaining their constant phase shift. Figure 2 
represents this oscillation for the free surface. 

As the effects of Coriolis acceleration were introduced 
into the model these conservative oscillations were greatly modified 
as can be seen from a comparison of Figures 2, 3* 11, and 12. The 
effect of a moving earth is that the total potential energy of the 
surface elevation will not be available for transfer to kinetic 
energy since the water particles are forced to turn due to the 
Coriolis Acceleration. Although the energy balance still remains 
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Figure 11. Effects of Inertial Rotation on Surface Elevation (PCD = 0) . 
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conservative, the relationship between the free surface elevation 
and the flow changes significantly. 

To understand how rotation acts to modify the flow, con- 
sider first the physical processes which drive the system. These 
are discussed in terms of the frictionless irrotational seiche. 

The initial conditions for the model were identical for all runs; 
that of a positive set-up of the free surface at the left boundary 
(x - 0) and no motion in the basin (a condition of static equili- 
brium). As the system is released the pressure gradients developed 
by this set-up start the fluid moving in the positive x-direction 
resulting in a net divergence in the left side of the channel and 
convergence in the right. The free surface responds to this flow 
by slumping at the left and rising at the right with equal magni- 
tudes. At time (t = T/4) the surface is level (i.e. no gradients) 
and the kinetic energy associated with the flow is equal to the 
potential energy of the initial free surface set-up. Therefore the 
surface continues to fall at the left and rise at the right until 
this kinetic energy has been entirely converted into potential 
energy and a state of momentary equilibrium is once again reached. 
The sea surface is now set up at the right boundary equal in magni- 
tude to the initial conditions, the pressure gradients equal but 
opposite to the initial values, and the flow reverses direction 
and the system returns to initial conditions, etc. It is important 
to realize that only the divergence of the flow in the cross-channel 
direction can influence the system because no gradients are allowed 
to develop along the length of the channel. 
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With these processes in mind the observed effects of ro- 
tation on the system can adequately be explained. With the initial 
set-up previously described, the flow will start in the same manner. 
As soon as the motion begins however, the Coriolis acceleration will 
start to turn the flow in the negative y-direction (left rotation), 
and a component of the flow will start to develop along the axis of 
the channel. Therefore some of the potential energy is converted 
into rotational kinetic energy. The amount of rotational kinetic 
energy that the system receives is a function of the velocity of 
the flow, the value of Coriolis, and the time scale of the motion. 
The slumping of the surface elevation at the left and rising at the 
right will continue as long as there is divergence in the positive 
x-direction; however since a certain amount of energy is rotational, 
the set-up at the right boundary can not achieve the magnitude of 
initial conditions. When the cross-channel flow vanishes the en- 
tire motion is in the along channel direction, with subsequent 
rotation turning the flow back in the direction from which it 
initially came. It should be noted that for a maximum value of 
(PCF = 1 ), Figure 3 shows the free surface elevation at the left 
boundary never falling below the horizontal. 

An interesting relationship developed in the transport as 
rotation was introduced into the system, resulting in a rectifi- 
cation of the along channel component of the transport (V) . An 
examination of Figure 12 showed that the cr oss -channel transport (U) 
periodically oscillates between maximum positive and negative values 
while the along channel component oscillated between zero and a 
maximum in one direction. It was easily demonstrated that this 
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direction is strictly a function of initial set-up. As the flow 
was initiated in the positive x-direction the initial along channel 
component developed in the negative y-direction. When the cross- 
channel flow reverses direction, along channel components are 
generated in the positive y-direction. These components are just 
sufficient to cancel the components previously established in the 
opposite direction, resulting in a rectification in the along- 
channel transport. For motions of large time scales, this process 
implies that important changes in the flow are caused by Coriolis 
accelerations and these will be discussed in the conclusions. 

Another observed feature of rotation was the subsequent 
reduction in oscillation period for increased values of Coriolis. 
The effect of rotation in reducing the spacial scale of the motion 
allows the system to oscillate with greater frequency. 

It is apparent from the previous discussions that the 
effects of inertial rotation can not be neglected when dealing with 
large time scale motion. 

3 . Combined Effects of Friction and Coriolis 

The remainder of the test runs were used to investigate 
how friction and inertial rotation in various combinations effect 
the transient motion. Some of the important consequences are de- 
picted in Figures 13 through 17. 

The phase shift difference between bottom stress and 
transport existed throughout all values of (PCF) and (PCD), as 
evidenced in Figures 13 and 14 . This phase lag occurred in varying 
degrees but was a definite characteristic of time dependent motion, 
vanishing only when the motion became critically damped. 
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Figure 13* Comparison of Bottom Stress and Mid-channel Transport (PCF = .5, PCD 
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Figure 14. Comparison of Bottom Stress and Mid-channel Transport (PCF = 1, PCD 
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Figure 15. Effects of Friction and Coriolis on the Free Surface 
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Figure 17. Hodograph of Mid-channel Velocity ( t - 5T) . 



The rectification of the along channel flow was also a pre- 



dominate characteristic of the flow where inertial rotation was 
considered. Figures 13 and l4 again show this relationship, as 
modified by friction, and Figure 15 depicts the combined effect on 
the free surface. It can be seen that for an extreme case (PCD = 1, 
PCF = 1) the free surface only very slowly slumps to a zero potential 
because friction retards the cross-channel flow long enough for ro- 
tation to reverse its direction before the free surface reaches the 
hor izontal . 

Figure 16 indicates that with rotation the energy level is 
damped with relatively fewer oscillations. It must be realized 
however that as (PCF) increases, the model is describing motion on 
an increasing time scale. 

Figure 17 represents the path that would be taken by fluid 
particles (or suspended matter) at three levels in the column over 
five periods of oscillation. In addition to the rectification pre- 
viously noted, there appears to be a net drift of the motion in the 
initial cross-channel transport direction. A closer examination of 
Figures 13 and 14 also show this drift and is a combined result of 
rotation and friction. Since the motion is being continuously 
damped by friction, the cr oss -channel transport can not return the 
fluid to its initial position after each oscillation and a net 
transport in the initial direction of motion results. The net 
drift is then strictly a function of initial set-up conditions of 
the free surface. The velocity profiles depicted in Figures l8 
through 20 show some of the fine details of the flow with the 
combined influences of friction and inertial rotation. 
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Figure 18. Effects of Friction and Inertial 
Rotation on Mid-channel Velocity 
(PCF = 1, PCD = 1) . 
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Figure 19. Effects of Friction and Inertial Rotation 
. . on Mid-channel Velocity (PCF = 1, PCD = 1, 
t = 3T/4). 
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Figure 20. Effects of Friction and Inertial Rotation on 

Mid-channel Velocity (PCF = 1, PCD =1, t = T) 
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Normalized Depth (-z/h) 



The controlling influence on the stability of the water 
column continues to be the friction, however the interactions of 
the fluid are now more complicated due to the introduction of 
inertial rotation. Under the influence of increasing frictional 
dominance towards the bottom, the left rotation induced in the flow 
by Coriolis is increasingly damped, resulting in the familiar Ekman 
spiral. During periods of the motion previously described as 
smooth, the spiral remains fairly constant as the column rotates 
to the left. However, during the transition periods (Figure 19) 
when frictional influence is reduced due to decreasing flow velocity, 
the spiral unwinds from the bottom and var ing amounts of twisting of 
the column results. This is again evidence of the transitional 
nature of the fluid during these periods. It should be noted that 
the velocity scale of Figure 19 is an order of magnitude smaller 
than that of Figure l8, indicating how small the velocities of the 
column are during transition. As the flow once again becomes uni- 
form, the spiral is quickly re-established under the stabilizing 
influence of bottom friction, as depicted in Figure 20. 

It is interesting to note from Figure l8 that during the 
transition periods the along channel component of velocity does in 
fact change direction for brief periods. Because these periods are 
of such short duration, this velocity reversal along the channel does 
not contribute to any significant transport, in this direction. 

C. DISCUSSION OF THE GEOPHYSICAL SIMULATIONS 

All of the previously discussed processes characteristic of 
time dependent motion, were observed during runs 26 through 30 * 
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However, since the magnitude of these processes are Coriolis and 
frictional dependent, they occurred on a much smaller scale than 
observed for the test cases. 

The model generated the time dependent variables produced by 
the seiche on a time scale equivalent to the resonance period of 
the basin. For these geophysical simulations, the periods were 
all too small to allow inertial rotation to have any significant 
effect on the motion. However, for shallow seas of large hori- 
zontal extent, values of (PCF) approaching .25 are entirely pos- 
sible and test runs 1 through 4 showed that inertial effects would 
be significant and must be considered. Additionally, motion whose 
time scale is driven by external forces, independent of the re- 
sonance basin, are likely to be significantly altered by rotational 
dependence. For example the long period waves propagating up and 
down the continental shelf probably experience all the rotational 
characteristics associated with the time dependent motion described 
in the previous section. The observed phenomenon of transport 
rectification might occur in these systems and suggests a mechanism 
for modifying bed load transfer across the shelf. This process is 
a function of initial conditions and appears to require a period of 
static equilibrium just prior to initial motion. A region where the 
seasonal migration of atmospheric pressure systems are predominantly 
in one direction, might satisfy these requirements, with a resulting 
net rectification in one direction. Although this process is highly 
speculative, it certainly appears possible and a modification to 
the model to facilitate its study is suggested in Section V. 

The frictional effects observed during runs 26 through 30 were 
also minimal. It should be noted that Ekman depths of greater than 
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20 meters are possible and thus these test runs represent a minimum 
fr ict ional dependence. 

The exception to these minimal effects was the observed phase 
lag between transport and bottom stress. This lag was pronounced 
in all cases and for the simulation of South San Francisco Bay, 
for example, averaged 40 degrees. This indicates that even small 
amounts of friction are important and re-emphasizes the need for a 
mathematical representation of time dependent transients that allows 
independent calculation of variables such as bottom stress and 
transport . 
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V. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER WORK 



The model demonstrated that eq. (22) was a suitable mathematical 
representation of the time dependent motion described by the governing 
equations and the initial/boundary conditions. In this formulation 
the time dependent variables were uniquely determined by the local 
time histories of the surface at each point and this allowed inde- 
pendent representation of the various transients. The model further 
demonstrated the feasibility of developing a totally implicit finite- 
difference scheme for approximating eqs . (3) and (22) with accept- 

able accuracy. With the insight gained from the results it should 
be possible to improve the accuracy of the system by removing some 
of the important restrictions imposed by the Ekman dynamics and 
the simplified model. 

The physical processes observed with time dependent motion 
showed the relative importance of friction and inertial rotation 
in determining the nature of the variables. It was apparent that 
for the larger time scales (PCF on the order of .25 or larger), 
Coriolis accelerations are important to the flow and must be con- 
sidered. However, shorter time scales will probably be adequately 
represented neglecting rotation with little accuracy loss. During 
the data analysis it appeared that friction, even in small amounts, 
exerted considerable influence on the characteristics of the vari- 
ables and their relationships, and should therefore not be neglected 
in representing time dependent motion. 
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Since frictional effects are predominant in the system, the 
simplified assumption of a constant coefficient of vertical friction 
throughout the layer, is probably one of the most limiting factors 
on the model’s accuracy. A technique for improving the repre- 
sentation of the influence of friction might be the assumption of 
constant stress layer (Prandl layer) underlying the Ekman layer. 

The solutions in these layers would be matched at their interface; 
that is, at some height above the bottom determined as a function 
of the amount of friction present and the roughness of the bottom 
boundary. This method has been used in air -sea boundary layer 
models with considerable improvement in accuracy, and merits con- 
sideration at the ocean boundary in this scheme. 

The assumption of constant depth does not account for modifi- 
cations of the flow caused by changes in the potential vorticity 
induced by the shrinking and stretching of the water column over 
varying topography. This restriction is probably most serious in 
shallow regions where time dependent flow is significant and changes 
in the bottom extreme, for example the flow in the vicinity of the 
Monterey Submarine Canyon. 

Certain large scale motions merit representation and investi- 
gation using this formulation. One of these is the study of 
seiches, especially the trapped modes, on the continental shelf 
since their buildup can cause storm surge (Pierson and Neuman, 

1966). Also, as previously mentioned, the large scale oscillations 
on the continental shelf might make important contributions to 
sediment transfer across the shelf. Therefore it might prove use- 
ful to remove the horizontal boundary conditions from the model 
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and drive the flow by an externally specified time dependent 
forcing term. This would indicate how effectively the motion was 
coupled to these terms and show whether the previously noted trans- 
port rectification is of sufficient magnitude to be significant. 

An examination of Tables I and II indicates how costly this 
type of numerical representation is in terms of computer time. It 
was not determined how much loss in resolution would occur by in- 
creasing the spacial grid size and reducing the length of the local 
time history considered. It is clear that considerable reduction 
in computing time would result and that such considerations are 
worth investigating before a similar model is used for extensive 
production runs. 
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APPENDIX A 



DETAILS OF THE FINITE- DIFFERENCE SCHEME 
Consider first the calculation of the components of the vector 



D . , recalling that 



d. = a©" + e 1 ? - a©" + AtA 

3 J + 1 3 3 - 1 



(Al) 



where 



a = 



At 

55 ? 



o o 

0 0 

-1 0 



C2/At 

D2/At 

0 



(A2) 



and 



e n = 

3 



A = 




{ d? 



(A3) 



(a4) 



D. = 

3 



(A5) 



Expanding eq. (Al) using (A2) through (A4) yields: 



72 



or , 
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eq. (49) and express it in inverse form as 



\ 



J 



(A6) 



(A7) 



(A8) 



-1 



^ = (I + o®.^) a, 3 = l 



(A9) 



is the (3 X 3) identity matrix. Using the appropriate 
for I, a, and ^ , it is easily showed that 
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(A10) 



(I + aE j _ 1 ) = 
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Therefore, 
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Where DET(y^ ^ ) is the determinant of the matrix (y^ and 
ADJ(y J is its adjoint. The computation of eq. (All) are ex- 
tremely lengthy and are not included. The resulting product of 
applying eq. (All) to eq. (A9) is as follows: 



E j DET_. 
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This completes the calculations of the components of E ^ . 
Finally consider the calculations of the components of the vector 
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F . = 
5 




Recall eq. (50) and express in its inverse form as 

F. = (I + aE. .. ) -1 (D.- aF. , ) , j S: 1 
D 1 J-l J J-l 



Again the computations of eq. (All) are involved and are not 
included. The results of applying eq. (All) to (A15) are as 
follows : 
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Note that the middle term of the vector is the entire quantity 
in the brackets. 

The implicit scheme for computing U, V, and § can now be 
constructed. Writing eq. (46) as 



^ n+1 +1 

e. . = E. n 9. + F . _ 

J - 1 3-1 3 D-l 



(A18) 



and substituting the quantities previously derived for E_. ^ and 
F 4 ^ • 

3-1 



yields : 
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Equations (A7) , (A12), (A13) , (Al6), (A17) , and (A19) through 
(A21) represent the expressions used to implicitly calculate values 
for at each grid point in the manner described in the text. 
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APPENDIX B 



DETAILS OP REPRESENTING Cl*, C2*, VELOCITY AND BOTTOM STRESS 



Consider the representation of Cl by first noting that the 
limits of integration have the following meaning: 
t = 0 represents the initial condition time, 
t = t represents the present time. 

t = t + At represents the time at which the value of the 

various variables are to be calculated. 

The integration of eq. (48) was performed numerically by 

dividing the total time into small intervals (At) and evaluating 

^(t ? ) as a constant at the mid-point of each interval. This 
* 

allowed Cl to be numerically represented as follows: 
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(Bl) 



With the further assumption that the effect of inertial ro- 
tation was constant in each interval, e ^ was approximated 

at the mid-point as e * 2 and treated as a constant. 

Analytical integration of the remaining integral resulted in 
eq. (51). To solve for Cl and D1 the following relation was used, 
and the real and imaginary terms collected, 
if t 



= cos(ft) - i sin(ft) 



(B2 ) 
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This allowed Cl and D1 to be represented as 
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Using the same assumption regarding constant inertial rotation 

-X- 

in each time interval, C2 readily reduced to eq. (52). Applying 
eq. (B2) to eq. (52) yields: 
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Consider now the scheme formulated to recover the velocity 
profile. Since velocity was calculated at time t + At, eq. (26) 



was represented by: 



w = 
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Applying the same assumptions and relations to eq. (B7) as 
were applied to eq. (Bl) yields: 
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wher e 



M 



A5 = (-1) cos(pz) 



a6 = (2£) 
v Ox 7 



The term (a 5) represents the depth dependence and (A 6) the 
slope in the most recent time interval. (A6) was calculated and 
stored by the model after U, V, and § were computed. 

The final quantity needing representation is bottom stress (T). 
Recalling eq. (54) and noting that the integral is identical to 
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that of eq. ( 26 ), bottom stress is easily represented as: 
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Wher e (T x ) and (Ty) are cross-channel and along-channel components 
of stress respectively. 

During the calculation of Cl, D1 , C2 , D2 , the model stores 
those common quantities needed for velocity and stress calculations 
as defined by eqs . (B3) through ( B6 ) , and (B8) thr ough (Bll ) . The 
program recalls these quantities and applies them to eqs. (B8) 

through (Bll) whenever it is required to calculate velocity and 

✓ 

bottom stress. 
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APPENDIX C 



COMPUTER PROGRAM DESCRIPTION 

The computer program was written in FORTRAN IV and was used 
with the IBM 360/67 computer system at the W.R. Church Computer 
Center, Naval Postgraduate School. The program has been generalized 
such that basin size, run time, initial free surface elevation, and 
number of modes considered in the solution can be varied in the data 
statements. Additionally, the model allows variable grid spacing, 
requiring only changes in the iraput parameters and DIMENSION state- 
ments. The influence of inertial rotation and vertical friction are 
introduced into the model by the two non dimensional ratios dis- 
cussed in the text. 

FORTRAN IV symbols used for the major model parameters are 
defined below. 

A Model time (real time normalized to the free oscillating 

seiche per iod) 

DH Spacial step in depth (fixed at one meter) 

DX Spacial step in acr oss -channel direction 

DT Time step 

ETA Free surface elevation 

ETAX Sea surface slope 

ETAXS Array for storing the time history of the sea surface 

slopes 

F Coriolis parameter 

FREQ Frequency of the free oscillating seiche 
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Acceleration of gravity 



G 



H 


Depth of water 


HI 


Free surface elevation at left boundary 


MAXH 


Integer representation of water depth (must be even) 


MAXL 


Maximum number of cross-channel grid points (must be odd) 


MAXP 


Interval in time step when velocity, energy, and bottom 
stress is calculated 


MAXT 


Maximum number of time steps 


PCD 


Ratio of Ekman depth to depth of the water 


PCF 


Ratio of Coriolis parameter to depth of the water 


pct 


Number of time steps per free wave period 


PE 


Potential energy of the seiche 


PERD 


Period of a free oscillating seiche 


R 


Coefficient of vertical friction 


RHO 


Density of the water 


SE 


Kinetic energy associated with the y-coraponent of velocity 


SIGMA 


Angular frequency of a free oscillating seiche 


ST R ESC 


y-component of bottom stress 


STRESR 


x-component of bottom stress 


T 


Real time 


TE 


Kinetic energy associated with the x-component of velocity 


TOT 


Total energy balance of the seiche 


UT 


x-component of volume transport 


VELC 


y-coraponent of velocity 


VELR 


x-component of velocity 


VT 


y-component of volume transport 
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The procedure for initialization of the program is as follows: 

1) Modify the initial DATA statements to correspond to the 
parameters of the motion being represented. 

2) Change DIMENSION statements accordingly. 

3) Insure that appropriate storage and computer time require- 
ments are satisfied. The program required a maximum of 100K BYTES 
storage on the IBM 360/67 system and variable CPU usage time as 
depicted in Tables I and II. 

4) A minor modification to the program is necessary when 
frictionless motions are being represented. If (PCD = 0) or 

(PCF = 0), remove the following three expressions from the program. 
K = ( I -1 ) - 10.0/CC 
IF(K.LE.0)K = 1 
DO 40 J = K, IMl 
Replace these expressions by 
DO 40 J = 1, IMl 

5) Because of the formulation, the program can not compute 
velocity at I = 1; therefore a value of MAXP = 1 should not be 
used. To recover the velocity profile at each time step (I ^ 1), 
initialize the DATA with MAXP = 2 and insert the statement 
(MAXP = 1) immediately before statement number 300 . 

The program, as appearing in the next section, is initialized 
in the MKS system and changes in the iraput should be made ac- 
cordingly. The major computational sections of the program are 
prefaced by appropriate comment statements, and a reader with 
FORTRAN IV experience should have little difficulty following the 
program . 
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COMPUTER PROGRAM 



DIMENSION ETA1< 11) ,ETA2<11) .ETAX(ll) ,ETAX1 (11) , 
1ETAXS(1 1,2 01 ), UT1 ( 11 ) ,UT2( 11 ) ,VT 1( 1 1 ) , VE LR ( 1 1 , 31 ) , 
2VELC (11,31) ,FUNR(11 ,25) ,FUNC(11 ,25 ) , P ( 25 ) , P2 (25 ) , 

3P3< 25) ,P4( 25), PK( 25) ,3< 31) ,C< 31) ,TIN< 11) ,PIN< 11) , 

4 POT (11) ,E11 <10 ) ,E21 <10 ) ,E31 ( 10 ), E13< 10), E23( 10) , 

5E33( 10) ,D1( 10) ,D2(10),D3(10) ,F1 < 10) , F2 <10 ) . F3 < 10 ) , 

6 D5T ( 10 ) , CR 1 ( 1 1 ) , CC 1 ( 1 1 ) , STRE SR ( 1 1 ) , STRE SC ( 1 1 ) 

C INPUT INITIAL DATA 

C 

DATA MAXT , MAXL , MAX H , MAX P , MAX M/ 200 , 1 1 , 30, 5,25/ 

DATA WIDTH, HI ,PC T, DH » RH 0/5000. ,1. ,40. ,1 . , 1027 ./ 

CAT A PIE,G, PCF, PCD/ 3. 14 15927, 9. 8, .5,0.0/ 

C 

C ECHO CHECK DATA 

WRITE(6, 310 ) 

WRITE (6, 320) MAXT , MAXL, VAXH, MAX P ,M AXM 
WRITE (6, 330 ) WI DTH, H I , DH ,RHO , PI E ,G , PC T , PCF , PCD 
C 

DX=WIDTH/FLOAT (MAXL-1 ) 

H= FLOAT ( MA XH ) 

PERD=(2 .0* WIDTH )/SQRT ( G*H) 

FREQ = 1. 0/PEP.D 
S IGMA=2 .0*P I E* FREQ 
F=PC F#S I GMA 

R=(PCD*H)**2*F/( 2. 0=*PI E**2 ) 

CT = PERD/PCT 
Q1=DT/(4.0*DX) 

Q 2 = 4. 0*DX 
C3=-2.C*G/H 

WRITE(6,350)H,WIDTH,DT,DX,P£RD,SIGMA,F,R 

C 

C COMPUTE AND STORE INITIAL CONDITIONS ON ETA ,E TAX, UT, VT 

C 

DO 10 L=1 , MAXL 
X=(L-1 )*DX 

ETA1 ( L) =HI*COS ( PI S*X /WIDTH) 

ETAX (L ) = (- 1. C*P IE^HI/WIDTH)*SIN( PI£*X/WIDTH) 

UT 1 ( L)=0.0 
VT1( L) =0. 0 
10 CONTINUE 
A=0. 0 

ETAX( MAXL) =0.0 
C 

C PRINT INITIAL CONDITIONS 

C 

WRI T5 ( 6, 360) A 

WRIT 5 (6, 370) (ETAl (L),L = 1,MAXL) 

WRI TE( 6,380) (ETA X( L) ,L=1 , MAXL) 

WRITS ( 6 . 400 ) ( UT1 ( L ) » L=1 , MA XL ) 

WRITE (6 ,405 ) (VT1 (L) , L=l, MAXL ) 

C CALCULATE AND STORE MODAL CONSTANTS 

C 

DO 20 M = 1 , MA X M 

P(M)=< 2.0*FL0AT ( M ) -1 .0 )*P I E/ < 2 .0*H ) 

P2( M)=P( M)**2 
PK ( M )=P2(M)*P 
P4( M)=EXP( -PK ( M ) *DT ) 

C 5 = P K ( M ) 

C4=C5*DT 

I F ( C4 . LT . 1 .E-4 ) GO TO 15 
P3( M) =( 1 • O-EXP (-C4 ) )/C5 
GO TO 20 

15 P3(M)=( DT*(1 .0-(DT*C5/2.0)*< 1 .0-C5*DT/ 3.0 ) ) ) 

20 CONTINUE 
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CALCULATE CR2 ANC CC2 

FDT = F*DT 
C1=C0S( . 5*FDT) 

C2=S IN ( • 5*FDT ) 

CR2=0 .0 

DO 30 M=1,MAXM 
CR2=CR2 + P3(M)/P2(M) 

30 CONTINUE 

CC2=-C3*C2*CR2 
CR2=C3*C 1*CR2 



BEGIN TIME LOOP 

DO 300 1=1 , MAXT 

T=I*DT 

A=T/PERD 

IFd.NE.DGO TO 35 
DO 31 L=1 «MAXL 
CR1 ( L ) =0 .0 
CCKU =0.0 
31 CONTINUE 

SET LEFT BOUNDRY VALUES FOR MATRICES 

Ell (1 )=0.0 
E21( 1) =0. 0 
E31 ( 1 )=0.0 
E13 ( 1 ) =0 .0 
E23( 1 ) = 0. 0 
E33 ( 1 ) = 1 .0/ .95 1 
F 1 ( 1 ) =C • 0 
F 2 ( 1)=0.0 
F3 ( 1 ) =0 .0 
DET( 1) =1.0 
GO TO 65 

(I.NE.l) CALCULATE CR1 AND CC1 

35 DO 60 L=1 » MAXL 
C6 = 0. 0 
C7=0 . 0 

DO 50 M=1 • MAXM 
FACR 1=0. 0 
FACR2=0 .0 
FACC1=0 .0 
FACC 2=0.0 
CC=PK(M)*DT 
I M1=I — 1 

K=(I-1)-10.0/CC 
IF(K.LE.0)K=1 
DO 40 J = K » I Ml 
C4=ETAXS (L » J ) 

r c = Kvp/_ r r : > s / T-1 - I \ l 

FACR1=FACR1+C4*CQS(FDT*( I+.5-J) )*C5 
FACR2=FACR2-C4*COS(FDT*< I-. 5-J) )*C5 
FACC1 = FACC1~C4 A S IN ( FOT * ( I + .5-J ) >*C5 
FACC 2= FACC 2+04*= SIN(FDT^(I-.5-J) )*C5 
40 CONTINUE 

FACR1 = FACR1*P4 (M) 

FACC 1=FACC 1*P4( M ) 

C5=P3 ( M ) 

FUNR(L,M)=FACR1*C5 
FUNC ( L ♦ M ) = FACC 1*C5 
C8=1.0/P2(M ) 

C6=C6+ ( FACR1+F ACR2 )*C5*C8 
C7=C7+(FACC1+FACC2)*C5*C8 
50 CONTINUE 

CR1 ( L) =C3*C6 
CC 1( L ) =C 3*C 7 
60 CONTINUE 
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CALCULATE VALUES OF MATRICIES 

65 MM1=MAXL-1 
DO 70 L=2 » MM1 

01(L )=<CR2/Q2)*(ETA1( L+l) -ETAl (L-l ) )+UTl (L )+CRl (L ) 
D2(L)=( CC2/Q2)*(ETA1(L+1)-5TA1(L-1 ) J+VTl ( L)+CC1( L) 

D3 ( L) =-01* (UT1 ( L + l ) -UT1 (L-l ) )+ETAl (L ) 

03=0 1/DE T ( L-l) 

Ell(L)=Q3*CR2*e33(L-l)/Q2 

E13( L) =03* (CR2 /DT-CR2*E13 ( L-l )/Q2 ) 

E21(L )=Q3*CC2^E3 3( L - 1 ) /0 2 
E23(L)=Q3*(CC2/DT-CC2*513(L-1)/Q2) 

E31( L) =-03* ( 1.0+CR2*E31<L-1)/Q2 ) 

E33 ( L ) = Q3*CP 2*E1 1( L-l) /Q2 

DET(L) =1.0-(1.0/02 )*(DT*E13(L )-CR2*E31 (L ) )+(CR2*DT/Q2 
1**2)*(E33(L )*E11(L )-E31 ( L) *E13( L ) ) 

Q4=D1(L)-CR2*F3(L-1 )/02 
05=D2( L)-CC2*F3 ( L-l )/02 
06=D3( L ) +DT*F 1( L-l )/Q2 

FI (L) = ( ( 1 .0 -DT + E13 (L-l )/ 02 )*Q4- ( CR2*E33 ( L- 1) /Q2 )*Q6 ) / 
1D£T(L-1) 

C4=(-CC2*E31 (L-l )/02 + CC2*DT*(E31(L-l)*E13(L-l)- 
1E3 3( L-l) -5 11 (L-l ) ) /Q2**2)*Q4 
C5=(1.0+(CR2*E31(L-1)-0T*E13(L-1) ) /Q2-CR 2* DT* ( E3 1 ( L-l ) 
1*E13 ( L-l ) -S3 3 (L-l)* Ell (L-l) )/Q2**2)*Q5 
C6=(CC 2*E33( L-l) /Q2)*Q6 
F2(L)=(C4+C5-C6)/D2T(L-1) 

F3(L)=( <DT*E11 (L-l ) /Q2 )*Q4+( 1 .0+CR2*E31 ( L-l )/Q2 )*Q6 )/ 
lDET(L-l) 

70 CONTINUE 

COMPUTE UT 2 » VT1. ETA2 

ETA2( MM1)=F3(MM1 )/(1.0-E33 ( MM1 )/ .951 ) 

ETA2(MAXL)=ETA2(MM1)/. 951 

UT2 ( MM1 )=E13 (MMl )*ETA2 (MAXL l+FKMMl) 

VT1( MM1 ) =E23 ( MM1 )*ETA2 (MAXD+F2 ( MM1 ) 

UT2 (MAXL )=0.0 
VT1 ( MAXL ) =0 .0 
MM3=MA XL-3 
DO 80 L= 1» MM 3 
M=MAXL-L 
C4=UT2( M) 

C5=ETA2(M) 

ETA2 ( M— 1 ) = E3 1 (M-l)*C4+£33( M-l )*C5+F3(M-1 ) 

UT2 ( M-l ) =iil 1 (M-l )*C4+?13(M-1)*C5+F1(M-1) 

VT1( M-l ) =521 (M-l ) +C4+E23 ( M-l ) *C5+F2 ( M-l ) 

80 CONTINUE 

ETA2 ( 1 ) = 5T A2 (2 ) / .951 
UT2( 1) =0.0 
VT 1 ( 1 ) =0 .0 

COMPUTE ETAX1 

DO 85 L=2 j MM1 

ETAXKL ) = ( ETA2( L + l )-ETA2(L-l) ) /(2.0*DX) 

85 CONTINUE 

ETAX1 ( 1 ) =0. 0 
ETAX1 ( MAXL )=0.0 

RESET ETAXS ARRAY AND RESET ETA1 , ETAX, AND UT1 
DO 90 L=1 » MAXL 

ETAXSt L,I)=(STAX1( L) + ETAX(L) )/2.0 
UT1 ( L ) = UT2 ( L ) 

ETAl ( L ) =ETA2 ( L ) 

ETAX( L ) =ETA XI ( L ) 

90 CONTINUE 

PRINT OUT UT1, VT 1 » ETAl, AND ETAX 



86 



WRI TE (6 ,410 ) A 

WRITE! 6,420) (UT1(L) ,L=1,MAXL) 

WRITE! 6 , 421 ) (VT1 (L ) , L = 1,MAXL ) 

WRITE (6, 42 5) (ETA1 (L) ,L=1,MAXL) 

WRITE ( 6,429) (ETAX( L) ,L = 1.MAXL) 

C 

I F ( MOO (I ,MAXP) • EQ. 0 ) GO TO 100 
GO TO 300 
C 

C I IS A MULTIPLE OF MAXP, CALCULATE VELOCITY 

C 

100 CO 115 L=1 , MAXL 
C6=ETAXS( L ,1 ) 

DO 110 M=1 , MAX M 

FUNR(L,M)=FUNR(L,M)+C1*P3!M)*C6 
FUNC! L,M)=FUNC( L,M)-C2*P3(M)*C6 
110 CONTINUE 
115 CONTINUE 

DO 140 L=1,MAXL 

MH1= MAX H + l 

DO 130 J =1 , MH1 

Z=(J-1)*DH 

SI=-1.0 

C4=0. 0 

C 5=0. 0 

DO 120 M=1 , MAXM 
C6 = P( M) 

C7=C0S( C6*Z ) /C 6 
C4 = C4+S I*FUNR (L , M ) *C7 
C5=C 5+SI*FUNC( L,M)*C7 
SI=-1.0*SI 
120 CONTINUE 

VELR(L, J)=-C3*C4 
VELC(L, J )=-C3*C5 
130 CONTINUE 

140 CONTINUE 
C 

C CALCULATE BOTTOM STRESS 

C 

C6=C3*R 

DC 142 L=1 ,MAXL 

C4=0.0 

C5=0.0 

DO 141 M=1 , MAX M 
C4=C4+FUNR(L,M) 

C5=C5+FUNC(L,M ) 

141 CONTINUE 

STRE SR ( L ) =C 4*C 6 
STR5SC!L ) = C5*C6 

142 CONTINUE 
C 

C CALCULATE ENERGY BALANCE USING SIMPSON'S RULE 

C 

145 DO 165 L=1,MAXL 
DO 150 J= 1 , MH1 
B(J)=VELR(L,J)*#2 
C( J)=VELC(L,J)**2 
150 CONTINUE 
EVSUM1=0.0 
EVSUM 2=0.0 
ODSU M1 = 0 .0 
0DSU^2=0.0 
EN SUM 1=B ( 1 ) + B( MH1) 

ENSUM2 = C (1 ) +C ( MH1 ) 

MH 2=MH 1- 1 
DO 155 J=2,MH2,2 
EVSUM1=EVSUM1+B( J) 

E V SUM2=E VSUM2+C ( J) 

155 CONTINUE 
MH2=MHl-2 
DO 160 J = 3 , MH2 , 2 
ODSU M1 = 0DSUM1 + B ( J ) 
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0DSUM2 = 0DSUM2+C( J) 

160 CONTINUE 

T IN( L ) = OH/ 3.0* (ENSUM 1+ 2.0*0 D SUM 1+4.0*EVSUM1) 

PIN( L)=DH/3. 0* (ENSUM2 + 2.0*0DSUM2+4.0*5VSUM2) 

POT (L )=ETA1(L)£*2 
165 CONTINUE 
E V SUM1 =0.0 
EVSUM2=0.0 
EVSU M3 =0 .0 
0DSUM1=0. 0 
0DSUM2=0.0 
0CSUM3 =0 .0 

ENSUM 1=P0T(1)+P0T( MAXL ) 

ENSUM2 = TIN(1)+TIN(MAXL ) 

ENSUM3 = PIN(1 ) + P I N ( MAXL ) 

DO 170 L=2 ♦ MM1 , 2 
EVSUM1=EVSUM1+P0T (L ) 

EVSUM2=EVSUM2 + TIN( L) 

EVSUM3=5VSUM3+P IN( L ) 

170 CONTINUE 
MM2=MA XL-2 
DO 175 L = 3 » MM2 » 2 
0DSUM1 =0DSUM1+P0T ( L) 

OD SUM2=0D SUM2+ TI N( L) 

0DSUM3=0DSUM3+P IN (L ) 

175 CONTINUE 

PE=. 5 : *RH0*G*(DX/3. 0)* (ENSUM1+ 2. 0*0DSUMl+4 .0*EVSUM1 ) 

TE = . 5*RH0* (CX/ 3.0 )*( ENSUM 2 + 2.0*00 SUM 2+4.0*E VSUM2 ) 

SE =. 5*RH0* (DX/3.0)* ( ENSUM3+2 .0*0DSUM3+4.0 *EVSUM3 ) 
TOT=PE+TE+ Sc 

PRINT VELOCITY 

WRITE (6 1 430 ) 

DO 180 J =1 , MH1 
JM1=J- 1 

WRITE (6, 440) JM1, (V SLR ( L , J ), L= 1 , MAXL ) 

180 CONTINUE 

WR IT E( 6 » 43 5 ) 

DO 190 J=1 ,MH1 
JM1= J- 1 

WR IT E( 6,440 )JM1, (VELC(L,J) ,L = 1,MAXL ) 

190 CONTINUE 

PRINT PE»TE»SE»TQT, STRE SR, AND STRESC 

WRITE ( 6,450) PE , Te , SE,TOT 
WRITE! 6, 45 5) (STRESRtL ) ,L = 1,MAXL ) 

WRITE (6, 460) (STRESC(L) ,L=1,MAXL > 

300 CONTINUE 

STOP 

310 FORMAT ( ' 1 ' , ' SC HO CHECK OF INPUT DATA',/) 

320 FORMAT ( *0' ,514) 

330 FORMAT! • O' ,951 1. 3) 

350 FORMAT (• 1* , 'SOLUTION TO SEICHE PROBLEM WITH CORIOLIS', 
1' FOLLOWS' , / , 1 X , ' DEPTH = ' , F5 .1 ,1X ,' METERS ' , / , IX, 

2 'WIDTH= ' »F 7. l,lX.'MeT5RS',/,lX.«DT=',F5.1,lX,'SEC'./, 
31 X, • DX= ' , F6 .1, IX, 'METERS •,/, IX, ' P5RD= ' , F 7. 1 , 1 X, 'SEC* , 
4/, IX,' SIGMA =' ,cll.3,/»lX,'F=' ,511.3,/ ,1X,' R=' ,£11 .3, // 
360 FORMAT!' • , 55X , • IN I TI AL CONDI TI ONS ',/, 56 X ,' TI M5 = ' , 
1F4.1 , IX, 'SECONDS' ) 

370 FORMAT!' • , 55X , ' ETA ' /( 11E11. 3) ) 

380 FORMAT!' • , 55X , ' ET AX '/ ( 11E11 .3 ) ) 

400 FORMAT!' ' ,55 X ,' TRANSPORT (UT )'/( 11E11 .3 ) ) 

405 FORMAT!' ', 55X .' TRANSPORT! VT) •/( 11E1 1 .3 ) ) 

410 FORMAT! 'O' ,55X, 'NORMALIZED T IM£= ' , F8 .2, / ) 

420 FORMAT!' • ,55 X ,' TRANSPORT (UT I •/< 11511.3 ) ) 

421 FORMAT! • ', 55X ,' TRANSPORT! VT »•/( HE 11. 3) ) 

425 FORMAT!' • ,55X , • ET A ' / ( 11E11 .3 ) ) 

429 FORMAT!' • , 55X , ' ETAX • / ( 1 IE 11 . 3 ) ) 
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430 FORMAT ( • 0 • , 5 5X , ' VELOC I TY ( VELR ) ' , / , 3X , ' Z ' , / ) 

435 FORMAT ( 'O' ,55X , ' VSLGC ITY (VELC I • , / » 3X , • Z • , / ) 

440 FORMAT ( • ' » 2X. I 2 ♦ ( HE 11. 3) ) 

450 FORMAT <'0','Pe='»2X,515.7,8X f ' TE= ' ,E 1 5. 7 , 8X , • SE 
1E15.7, 8X T ' TOT=* ,E15.7) 

455 FORMAT ( . 55X. • STRESS(R) • ,/( 11E11.3) ) 

460 FORMAT ('O' , 55 X , • ST R ESS ( C ) • , / ( 1 IE 1 1 .3 ) ) 

END 
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